## Abstract

Rayleigh-scattering radiance (*L _{r}*) calculations based on the standard algorithm are often associated with significant uncertainties leading to inconsistent water-leaving radiance retrievals, both spatially and temporally across latitudes and altitudes. The uncertainty could result from the use of Rayleigh lookup tables generated for the standard surface atmospheric pressure and hence the Rayleigh optical thickness (ROT) at the specific atmospheric pressure regardless of its daily and seasonal variations. This study presents a new algorithm (hereafter referred to as the refined algorithm) to compute the Rayleigh-scattering radiance that relies on accurate calculations of the ROT as a function of the composition of air (CO

_{2}volume concentration), surface atmospheric pressure and relative air mass for given sun-sensor geometries. As CO

_{2}is well mixed throughout the atmospheric column, the CO

_{2}volume concentrations derived from this study agree well with measurements in different seasons across studied latitudes. Relative air mass has a significant effect on the ROT and that is calculated as a function of apparent sun-sensor zenith angles with the variations in pressure and thermal characteristics of the atmosphere. Thus, the results indicate significant variations of ROT and air mass with location on the earth’s surface and their influence on the

*L*, particularly in the UV-Blue region of the spectrum. The refined algorithm for calculating the

_{r}*L*is tested on several MODIS-Aqua Level 1A data and the relative errors in Rayleigh-scattering radiance and normalized water-leaving radiance (

_{r}*Lwn*) retrievals between the refined algorithm and standard (SeaDAS) algorithm are compared using in-situ measurement data collected at MOBY (clear ocean), AERONET (turbid coastal ocean), and NOMAD (clear ocean) sites. The results indicate that the

*L*calculated using the SeaDAS algorithm are mostly underestimated and show significant departures with the

_{r}*L*calculated using the refined algorithm. This departure induced by the SeaDAS algorithm to

_{r}*L*becomes larger with decreasing wavelength (Δ

_{r}*L*from −2.38% at 412 nm to 1.69% at 678 nm), which causes errors in

_{r}*Lwn*retrievals (Δ

*Lwn*) of up to 26.48% at 412 nm and 13.34% at 678 nm. The overall improvements in the retrieved

*Lwn*values achieved vary from 56% at 412 nm to 29% at 678nm, which yield similar improvements in

*Lwn*retrievals with lower errors and higher slopes and correlation coefficients when compared with the in-situ

*Lwn*data. These results indicate that the refined algorithm for computation of the

*L*can yield more accurate

_{r}*Lwn*retrievals and produce spatially and temporally consistent biogeochemical products at different latitudes and altitudes as desired by the scientific community.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Water color observations provided by the space borne sensors – such as SeaWiFS (Sea-Viewing Wide Field-of-View Sensor), MODIS-Aqua/Terra (Moderate-Resolution Imaging Spectroradiometer), MERIS (Medium-spectral Resolution Imaging Spectrometer), OCM-1 and 2 (Ocean Color Monitor), GOCI (Geostationary Ocean Color Imager), VIIRS (Visible Infrared Imager Radiometer Suite), Hyperion, HICO (Hyperspectral Imager for Coastal Oceans), Sentinel-3 (Ocean Land Color Instruments), and OLI (Operational Land Imager) – have become an exceptionally powerful tool for quantitative assessments and monitoring of coastal and inland water ecosystems and achieving the expanding needs of the research community for applications including biogeochemical variability, carbon cycle, water quality, primary productivity, phytoplankton functional types, harmful algal blooms (HABs), river blooms, sediment dynamics, oil spills, optical bathymetry, bottom detection and classification, and coral reef and sea grass mapping [1–14]. These sensors measure the top-of-atmosphere (TOA) signal from the surface-atmosphere system at a number of discrete bands (typically 10-40 nm) or narrow bands (typically 10 nm or better) in visible and near-infrared parts of the spectrum. The major portion of the TOA signal (∼80-90%) is contributed by the atmospheric (path) scattering signal due to air molecules and aerosol particles. The process of removing the atmospheric path signal from the TOA signal in order to retrieve the desired water-leaving radiance signal is referred to as atmospheric correction [15]. The most dominant signal in the surface-atmosphere system for visible bands is from the contribution of scattering by air molecules (i.e., Rayleigh-scattering radiance, *L _{r}*), and this signal contribution progressively increases toward shorter wavelengths because of the inverse relationship of the scattered intensity to the wavelength of radiation. Notably, molecular scattering is the strongest component of the atmospheric extinction of solar radiation for visible bands in the cloudless atmospheric conditions over the oceanic and land areas (including Arctic and Antarctic regions) [16–18]. This scattering component is often dominated by gas molecules such as O

_{2}, N

_{2}, CO

_{2}, and Ar [18]

*.*Because the Rayleigh-scattering radiance is the most dominant component of the TOA signal, accurate models are required to calculate the Rayleigh-scattering radiance contributions in visible wavelengths in order to retrieve consistent water-leaving radiance and biogeochemical products spatially and temporally.

In the atmospheric correction procedure developed for the Coastal Zone Color Scanner (CZCS), Gordon et al. [19] developed an algorithm to estimate the *L _{r}* for visible bands assuming a flat ocean surface with the inputs of solar-sensor geometries. For a wind-roughened ocean surface, [20] improved the algorithm for atmospheric correction of ocean color sensors that estimates the

*L*as a function of the ocean surface wind speed (which causes sea surface roughness) and specific atmospheric pressure for given solar-sensor geometries. These calculations are based on the Rayleigh lookup tables generated for the standard atmospheric pressure and the variation of the atmospheric pressure over the water bodies at different latitudes and altitudes (high altitude lakes) is compromised. To overcome this limitation, further modifications were applied to the previous algorithm with variation of the atmospheric pressure at surface level [17]. Although the refined pressure correction algorithm showed noticeable improvements in the satellite-retrieved water-leaving radiances for visible wavelengths, uncertainties still remain in the

_{r}*L*computations over the oceanic and inland water systems at different latitudes and altitudes. The earlier calculations of Rayleigh-scattering optical thickness (ROT) using the atmospheric pressure data showed a variation of up to 1% in ROT values [21]), which could have an effect on the ROT computations as a function of the standard atmospheric pressure and hence on the

_{r}*L*calculations for visible bands.

_{r}Rayleigh-scattering optical thickness values obtained from earlier studies [22] were based on a fitting equation that computes the ROT as a function of wavelength for a specific atmospheric condition. However, recent studies showed that Rayleigh scattering varies as a function of wavelength with respect to the atmospheric conditions (refractive index changes as a function of wavelength and altitude depending on air pressure, air temperature, and partial pressures of the atmospheric constituents N_{2}, O_{2}, Ar, CO_{2} and H_{2}O) [18,23]. On the other hand, calculations of the values of relative air mass using the secant approximation (as a function of apparent solar zenith angle) often contain the uncertainties that could affect the total ROT values at higher latitudes under low solar elevation angles.

Considering the accuracy required for Rayleigh-scattering radiance calculations, the present study is aimed to develop a refined algorithm to calculate the *L _{r}* at UV, visible and IR wavelengths based on improved parameterizations accounting for the variations in the physical and compositional effects of the atmospheric system. Accordingly, the Rayleigh-scattering radiance calculations depend on the surface pressure, atmospheric constituents, solar-sensor geometries, surface wind speed, and sea surface roughness. Unlike the SeaDAS algorithm, the ROT is calculated as a function of the surface pressure, CO

_{2}volume concentration, and relative air mass for given solar-sensor geometries. The Rayleigh-scattering radiances calculated using the regenerated Rayleigh lookup tables are then compared with those derived using the SeaDAS algorithm for a variety of the oceanic and atmospheric conditions recorded by the MODIS-Aqua sensor. For the accuracy assessment, concurrent MODIS-Aqua and in-situ measurements were obtained over the MOBY, AERONET-OC [24] and NOMAD sampling locations (Fig. 1, description in a later section), which cover a wide range of oceanic conditions (clear open ocean to turbid coastal waters) and atmospheric conditions (low to high latitudes). To evaluate the accuracy of the

*L*and

_{r}*Lwn*products derived from satellite data using these algorithms, some important statistical matrices are used to quantify the deviations and errors in the satellite derived products.

## 2. Theoretical background

The satellite sensors measure the spectrum of reflected light radiance at the top-of-atmosphere (TOA) which can be written as

*L*and

_{TOA}(λ), L_{r}(λ), L_{a}(λ), L_{ra}(λ), L_{g}(λ), L_{wc}(λ)*L*are the TOA radiance, Rayleigh-scattering radiance, aerosol-scattering radiance, Rayleigh and aerosol-scattering radiance (interaction term), sun glint radiance, whitecap radiance, and water-leaving radiance, respectively [17,25].

_{w}(λ)*T(λ)*and

*t(λ)*are, respectively, the direct and diffuse transmittances along the sensor viewing direction. In the atmospheric correction procedure, the contributions of the atmosphere and sea surface (specular) reflection are removed from the TOA signal in order to retrieve the desired water-leaving radiance. To avoid miscalculations of

*L*and

_{g}(λ)*L*under high wind speed conditions, these contributions are generally ignored for brevity. The contributions

_{wc}(λ)*L*and

_{a}(λ)*L*are due to aerosols and multiple scattering effects and can be estimated using the Gordon and Wang [25] algorithm when combined with the NIR iterative scheme [26] or NIR-SWIR scheme [27]. The aerosol radiances estimated for the NIR or SWIR bands are then extrapolated over visible bands based on the selected aerosol models (Ahmad et al., 2010). For optically complex oceanic and inland waters, a Gaussian extrapolation method using the spectral shape parameters [1] has been successfully demonstrated using MODIS-Aqua, Landsat-OLI and HICO images [28–30].

_{ra}(λ)Because *L _{r}(λ)* is the most dominant signal in the

*L*, it is critical that the Rayleigh-scattering radiance must be computed accurately. According to Gordon et al. [19],

_{TOA}(λ)*L*was computed using the single scattering approximation for analysis of CZCS imagery,

_{r}(λ)*F*= instantaneous extraterrestrial solar irradiance,

^{’}_{0}*τ*= Rayleigh-scattering optical thickness,

_{r}*ρ(θ*= Fresnel reflectance of the interface for an incident angle

_{0})*θ*, and

_{0}*P*= Rayleigh scattering phase function.

_{r}(α, λ)*θ, ϑ, θ*, and

_{0}*ϑ*, are the zenith and azimuth angles of a vector from the sea surface to the sun and sensor, respectively [19].

_{0}*P*(

_{r}*α*) is given by,

*I*is the extraterrestrial solar radiation flux at wavelength

_{0}(λ)*λ*,

*I(λ)*is the flux reaching the surface, and

*τ(λ)*is the total optical depth. The total optical depth of the atmosphere is defined by several components, ascribed to Rayleigh scattering, aerosol scattering, and gaseous absorption [31]. Considering Rayleigh-scattering radiance calculations, estimates of ROT are often based on the fitting equation using Penndorf data [32]. Earlier studies computed ROT values accurately over a particular wavelength range (e.g., visible range) and results of such calculations are inaccurate over the other wavelength range (e.g., UV range) due to the limitation of curve fitting techniques [33]. For these wavelength ranges, the conventional power-law equation used to derive the ROT is significantly different because of the constant slope (‘

*B*’ in Eq. (12) of Bodhaine et al. [31]) regardless of its variation in the visible-NIR region. Some studies have suggested different formulations to calculate ROT values over a greater range of the EM spectrum [34]. Hansen & Travis [22] suggested an equation of the form

*τ*is the ROT using the standard atmosphere pressure. For the atmospheric pressure variation at a site, ROT is defined as where

_{o}(λ)*P*= site-specific atmospheric pressure, and

*P*= standard atmospheric pressure (normalized to 1013.25 mb). In the atmospheric correction procedure, Rayleigh lookup tables are generated for the standard atmospheric pressure as a function of wavelength, the solar–sensor geometry, wind speed (sea surface roughness), and atmospheric pressure [19,35]. According to Gordon et al. [19], the Rayleigh-scattering radiance changes due to the variation of the sea surface atmospheric pressure can be determined using

_{0}*θ*is the zenith angle,

*τ*and

_{r}(λ, P_{0}+∇P)*τ*are the Rayleigh optical thicknesses, at a given wavelength, for atmospheric pressures of

_{r}(λ, P)*P*and

_{0}+∇P*P*respectively. Later, Wang [17] redefined the Gordon [19] algorithm as

_{0}*M = 1/cos θ*is the air mass value for the solar and sensor zenith angles of

_{0}+1/cos θ*θ*and

_{0}*θ*, and

*C (λ, M)*is the attenuation coefficient from the fitting equation. This solution provides accurate calculations of Rayleigh-scattering radiance at a given wavelength as a function of air mass and solar-sensor geometries for standard atmospheric conditions. Because the ROT varies spatially, temporally and vertically depending on the physical, thermal and compositional characteristics of the atmosphere, this scheme provides inaccurate and less consistent results at different latitudes and altitudes (high altitude lakes).

## 3. Description of the refined algorithm

In this procedure, the Rayleigh-scattering radiance (*L _{r}*) is computed for a given wavelength by accounting for the variations of the CO

_{2}volume concentration, pressure and air mass within the atmospheric column. Figure 2 describes the implementation of this scheme for

*L*computations. The refined formula takes the form

_{r}*P*is the pressure at a site,

*P*is the standard pressure (1013.25 mb), and

_{0}*L*is the Rayleigh-scattering radiance for the standard atmosphere.

_{r}(τ_{0}(λ, P_{0}))*C*is the dynamic constant derived from the air mass and optical thickness of the atmosphere as where

*τ*is the ROT calculated as a function of wavelength and CO

_{r}_{2}volume concentration, and

*M*is the relative air mass for the solar-sensor zenith angles. To facilitate this calculation, Rayleigh-scattering optical thickness must be determined for a given wavelength and composition as it is dependent on the atmospheric pressure and composition at the site. Based on the

*τ*calculated for a given wavelength and composition, the Rayleigh lookup tables can be regenerated with the new values of ROT. Following Penndorf [32] and Bodhaine et al. [31], the Rayleigh scattering coefficient is calculated from the equation

_{r}*σ*is the scattering cross section per molecule (cm

^{2}molecule

^{−1}),

*N*is molecular number density per unit volume,

_{s}*n*is the refractive index of the gas, (6 + 3ρ)/ (6 -7ρ) is the depolarization term or F(air) or the so-called King factor. The King factor depends on the gas mixture,

*N*depends on temperature and pressure, and the resulting

_{s}*σ*depends on the composition of the gas but is independent of temperature and pressure [31].

The refractive index of standard air at 300 ppm of CO_{2} concentration is calculated accurately using the equation of Edlén [36]:

_{2}is given by The depolarization of air is estimated as a function of the atmospheric constituents such as O

_{2}, N

_{2}, Ar and CO

_{2}[37] as

_{2}(C

_{CO2}) is expressed in parts per volume by percent, and x is a correction factor for other gas concentrations (according to CO

_{2}concentration). The depolarization of N

_{2}, O

_{2}, Ar and CO

_{2}are a function of the wavelength. According to Bodhaine et al. [31], the Rayleigh-scattering optical depth is determined for a given wavelength and composition (i.e., molecular number density per unit area in the atmospheric column at the site). The molecular number density depends only on the pressure, and hence, the ROT is calculated as where

*P*= pressure,

*A*= Avogadro’s number,

*m*= mean molecular weight of the air, and

_{a}*g*= acceleration of gravity. According to List [38], the sea level acceleration of gravity is based on the pixel location (latitude) and the mean molecular number density at the site. This can be calculated based on the CO

_{2}volume concentration as it is well mixed within the atmosphere column rather than the other atmospheric constituents.

Seasonal and latitude variations in atmospheric pressure and increasing trend of CO_{2} concentration in atmosphere are considered for these calculations. To account for these changes, it was decided to include the variability in refractive index, molecular number density, and King factor based on the CO_{2} concentration. Latitude-wise monthly mean atmospheric CO_{2} records of the National Oceanic and Atmospheric Administration (NOAA) Greenhouse Gas Marine Boundary Layer Reference (2000-2005) data have been used for formulating fitting equations for the CO_{2} concentration. The monthly average increment of CO_{2} concentration is calculated from the data for the period of 2000 to 2005 (as shown in Table 1). Figure 3 shows the comparison of NOAA monthly mean atmospheric CO_{2} concentrations and the corresponding atmospheric CO_{2} concentrations calculated from the fitting equations.

The results of the calculations clearly showed the strong seasonal variation in CO_{2} with latitude and the fitting equation estimated CO_{2} are in good agreement with monthly mean CO_{2} measurements. Based on the CO_{2} concentration, the refractive index of air, depolarization factor, and mean molecular number density are calculated using Eqs. (13) and (14). Then, the Rayleigh-scattering cross section is estimated using Eq. (11) and the air mass is estimated between the sun to surface and the surface to sensor using Eq. (16) [39]. Considering the high precision required by ROT calculations, relative air mass (*M(θ)*) is calculated as a function of the solar-sensor zenith angles and pressure variation at each stratified layer. Note that the vertical profiles of atmospheric molecular number density are influenced by the variations in pressure and temperature conditions and the atmospheric constituents [39].

According to Thomason et al. [40], the relative optical air mass of the attenuating atmospheric constituent is calculated as a function of apparent solar zenith angle *θ* as,

*P*is the standard atmospheric pressure, and

_{0}*n*and

_{0}*n*are the refractive indices of the lower and upper atmospheric layers, respectively.

_{z}*z*and

_{0}*z*are altitudes of the respective atmospheric layers, and

*∇P*is the pressure difference between the layer. Figure 4 illustrates the relative optical air mass variation (at 412 nm) as a function of the solar-sensor zenith angles above the horizon. Consistent with the earlier results [39,41], the air mass reaches close to 1 at sea level, when the sun is directly overhead (zenith angle = 0) (Fig. 4). As the solar zenith angle becomes larger, the air mass increases with the optical path length of the atmosphere. Similarly, as the sensor zenith angle becomes larger, the optical path length of the atmosphere grows larger and the air mass is increased. Using all of the above values of the atmospheric parameters, the calculations of τ

_{r}and

*L*are performed as accurately as possible because a few percent error in

_{r}*L*can translate into several percent error in the retrieved water-leaving radiance products.

_{r}## 4. Data sets

#### 4.1 In-situ data

The efficiency of the refined Rayleigh-scattering radiance algorithm is assessed by comparing satellite-retrieved products with in-situ measurements from the open ocean and coastal regions under various atmospheric conditions. The in-situ radiometric measurements made at the MOBY (Marine Optical BuoY at open ocean site) (https://www.star.nesdis.noaa.gov) and six AERONET-OC (Aerosol Robotic Network - Ocean Color) sites (https://aeronet.gsfc.nasa.gov) were used for this analysis. These six AERONET-OC sites were chosen as representative cases for coastal oceanic waters, including Gustav_Dalen_Tower (58.594N, 17.467E), Helsinki_Lighthouse (59.949N, 24.926E), WaveCIS_Site_CSI_6 (28.867N, 90.483W), Thornton_C-power (51.5332N, 2.955E), MVCO (41.325N, 70.567W), and LISCO (40.955N, 73.342W). Using these data, the normalized-water leaving radiance was determined for visible and NIR wavelengths $Lwn(\lambda ) = {F_0}(\lambda ) \times {{{L_w}(\lambda )} \mathord{\left/ {\vphantom {{{L_w}(\lambda )} {E{d_s}(\lambda )}}} \right.} {E{d_s}(\lambda )}}$, where *F _{0}(λ)*– Extra-terrestrial solar irradiance;

*L*– Water leaving radiance;

_{w}(λ)*Ed*– Downwelling solar irradiance.

_{s}(λ)For open oceanic waters across the latitudes (Chl < 0.5 mg m^{−3}), several in-situ measurements from the NASA bio-Optical Marine Algorithm Dataset (NOMAD) were considered for validating the algorithm results. The NOMAD dataset is a global, high-quality, in-situ, bio-optical dataset, publicly available for algorithm development and satellite ocean color data validation. This dataset comprises laboratory-analyzed biological parameters such as pigment concentrations including chlorophyll-*a*, chlorophyll-*b*, alloxanthin, peridinin and zeaxanthin, and in-situ measured radiometric quantities including spectral downwelling irradiance and water-leaving radiance data. In addition, latitude-wise monthly mean CO_{2} concentration of the atmosphere for the period of 2000-2005 were obtained from the National Oceanic and Atmospheric Administration (NOAA) Greenhouse Gas Marine Boundary Layer Reference database (https://www.esrl.noaa.gov) for deriving parameterizations for computation of Rayleigh-scattering optical thickness. The in-situ measurement data described above are sufficiently representative of a variety of the oceanic and atmospheric conditions.

#### 4.2 Satellite data

The refined Rayleigh-scattering radiance algorithm is tested on several MODIS-Aqua Level-1A data (with 1 km spatial resolution at nadir) acquired over the MOBY, AERONET-OC and NOMAD sampling locations. These MODIS-Aqua data were downloaded from the NASA Goddard Space Flight Center (http://oceancolor.gsfc.nasa.gov/) and processed to the Rayleigh-scattering radiance and normalized water-leaving radiance using the SeaDAS system (with refined algorithm). The effects of the refined algorithm on these satellite-derived products are assessed by comparison with those of the standard algorithm (SeaDAS). For retrieving the MODIS-Aqua normalized water-leaving radiance (*Lwn*), the Rayleigh-corrected radiances computed by refined and standard algorithms were input to a common aerosol correction algorithm [1,25] to output the *Lwn* products for visible and NIR wavelengths. Comparison of these results at different data product levels enables the importance of the refined algorithm to be put in context.

## 5. Results and discussion

The refined algorithm is based on the same physics as the SeaDAS algorithm, which is capable of yielding accurate Rayleigh-scattering radiances with the inputs of new ROT values. The ROT values computed using the refined algorithm are sufficiently accurate to accommodate the variations in atmospheric conditions induced by vertical and lateral variations in pressure, composition, and air mass over different locations and seasons, and thereby yield more consistent water-leaving radiance retrievals across the latitudes and altitudes. The refined algorithm was incorporated into the SeaWiFS processing system and applied to several MODIS-Aqua images acquired over the clear oceanic waters (MOBY and NOMAD) and turbid coastal oceanic waters (AERONET-OC). The results were systematically analysed by comparing with those derived from the SeaDAS algorithm and in-situ measurements. In addition, several test cases using MODIS-Aqua images from different locations and seasons were carried out to assess the relative differences of the refined and standard algorithms.

#### 5.1 Accuracy comparisons using measurement data

To see the improvements in the retrieved *Lwn* values due to the refined algorithm, the MODIS-Aqua derived *L _{r}* and

*Lwn*obtained using the refined and SeaDAS algorithms were compared with in-situ measurements at the MOBY, NOMAD and AERONET-OC sites. Figure 5 provides examples of the MODIS-Aqua derived

*Lwn*spectra for cases of clear atmospheres (i.e., more Rayleigh-scattering contributions) in clear open oceanic waters (MOBY) and moderately turbid, turbid, and phytoplankton-dominated waters (AERONET-OC – Thornton_C-power, WaveCIS, and Helsinki). When the aerosol correction scheme remained fixed in our atmospheric correction (SS14) scheme [2], it is seen that when computed using the SeaDAS Rayleigh-scattering radiances,

*Lwn*values were overestimated showing significant departures with the measured

*Lwn*values due to the underestimated

*L*values across visible wavelengths. For comparison purposes, the INIR scheme (available in the SeaDAS system which is commonly used to estimate aerosol radiances in global oceanic waters) was used to retrieve

_{r}*Lwn*values for clear oceanic waters at the MOBY station (results produced by the SeaDAS-INIR algorithm for the AERONET-OC data are not included due to the biased

*Lwn*caused by significant amounts of particulate matter). The results suggest that the SeaDAS-INIR scheme produced similar

*Lwn*values (as the Refined-SS14 algorithm) in clear oceanic waters (MOBY), which confirms the accuracy of aerosol estimation by the SS14 algorithm. However, validation of the SeaDAS-SS14 shows that the departure between the retrieved and measured

*Lwn*values in turbid and productive waters (AERONET-OC data) increases with decreasing wavelength by several times. Because phytoplankton pigment retrievals rely on spectral ratios in the blue-green bands [9], the poor (underestimation) approximation of

*L*computations by the SeaDAS algorithm is expected to seriously degrade (lower) the pigment retrievals and introduce inconsistencies in the derived spatial and temporal water color products. This problem will become more serious when applied to new generation sensors having UV bands because of the increasing difference between the retrieved and measured

_{r}*Lwn*values at the short wavelengths. In contrast, the more accurate computations of

*L*by the refined algorithm significantly reduced the discrepancy and produced good agreement between MODIS-Aqua and in-situ measurements of

_{r}*Lwn*under various atmospheric and oceanic conditions.

Statistical evaluation of the data based on the relative errors in Rayleigh-scattering radiance (Δ*L _{r}* =

*L*

_{r}_{SeaDAS}-

*L*

_{r}_{Refined}/

*L*

_{r}_{Refined}) and normalized water-leaving radiance (Δ

*Lwn*=

*Lwn*

_{SeaDAS}-

*Lwn*

_{Refined}/

*Lwn*

_{Refined}) and the mean relative error in normalized water-leaving radiance (MRE =

*Lwn*

_{Satellite}-

*Lwn*

_{in-situ}/

*Lwn*

_{in-situ}) gave further insights into the disparities between these two algorithms. Because there was the strong variation of the relative difference in the values of

*L*for visible bands between the SeaDAS and refined algorithms, the error Δ

_{r}*L*reached as much as −2.39% for the wavelength 412 nm and decreased to −1.69% for the wavelength 678 nm. In this case a few percent error in

_{r}*L*eventually translated into several percent error in the derived

_{r}*Lwn*. Considering these twelve samples (Fig. 5), the error Δ

*L*for wavelengths at 412, 443, 488, 555, and 678 nm were, respectively, −2.39%, −2.18%, −1.99%, −1.83% and −1.69%, which translated into the error Δ

_{r}*Lwn*64.42%, 41.48%, 19.11%, 9.96%, and 12.02% in the

*Lwn*products. Consequently, the SeaDAS algorithm yielded MRE of 89.4%,53.6%, 26.4%, 37.8%, and 159% at 412, 443, 488, 555, and 678 nm, respectively (higher MRE values at longer wavelengths due to the subsequent aerosol correction). Because the refined algorithm produced closely consistent

*Lwn*retrievals with in-situ measurement data, the MRE is significantly reduced to 51.8%, 31.1%, 19.3%, 31.1%, and 133.6% for the wavelengths mentioned above. Significant errors associated with computation of the

*L*(SeaDAS algorithm) seriously compromise the accuracy of atmospheric correction and degrade the quality of derived biogeochemical products from satellite data.

_{r}The accuracy of the refined algorithm was further determined based on more matchup data of MODIS-Aqua and in-situ measurements from the MOBY and AERONET-OC sites for cases of clear atmospheres. Figure 6 shows the MRE distributions across the samples wherein the Δ*L _{r}* calculated between the SeaDAS and refined algorithms were −2.38%, - 2.17%, −1.98%, −1.82%, and −1.69% at 412, 443, 488, 555, and 678, respectively (N = 56) due to the underestimation of

*L*by the SeaDAS algorithm. In this case, the error Δ

_{r}*Lwn*, corresponding to the error Δ

*L*, were 26.48%, 19.20%, 12.20%, 8.51%, and 8.66% at 412, 443, 488, 555, and 678, respectively. Note that the resulting sample-averaged MRE of the SeaDAS algorithm for these wavelengths were significantly higher (87%, 51%, 20%, 24%, 147%) than those (31%, 20%, 7%, 13%, and 118%) of the refined algorithm. Clearly, the error with the SeaDAS algorithm becomes progressively larger in the shorter wavelengths due to the underestimated

_{r}*L*values (Fig. 6).

_{r}The significant error in the computations of the *L _{r}* was then propagated through the atmospheric correction process and caused biased (overestimated)

*Lwn*values. In contrast, the refined algorithm achieved significant improvements in the retrieved values of

*Lwn*for the MODIS-Aqua wavelengths (i.e., MRE 56%, 31%, 13%, 11%, and 29% at 412, 443, 488, 555 and 678 respectively). Note that unlike the error Δ

*Lwn*, a fraction of the error MRE in the

*Lwn*retrievals could be caused by the differences in MODIS-Aqua pixel and in-situ point measurements due to sub-pixel and dynamical variability of the water constituents. Nevertheless, these matchup data exhibited lower MRE, RMSE and intercept values and higher slopes and correlation coefficients for the refined algorithm than for the SeaDAS algorithm (Table 2), which confirms significant improvements of the refined algorithm over the SeaDAS algorithm.

Figure 7 shows histograms accurately representing the frequency distributions of MRE in the retrieved values of *Lwn* using the refined and SeaDAS algorithms for the MODIS-Aqua and MOBY/AERONET-OC matchups in open ocean and coastal waters. The histogram representation of data shows that significant departures between the values of MRE using the refined and SeaDAS algorithms are seen progressively increasing with decreasing wavelength. It can be seen that most of the data are centered around the origin with a normal, symmetrical peak with fewer larger MREs (still lower) due to the refined algorithm. Because of significant errors induced by the SeaDAS algorithm, many of these data are slightly shifted to the right with a broadened offset peak indicating larger MRE values particularly at the blue wavelengths. This indicates higher performance of the refined algorithm and lower performance of the SeaDAS algorithm at visible wavelengths.

The algorithm performance was further assessed using MODIS-Aqua and NOMAD in-situ matchup data (several samples considered across various latitudes) from clear oceanic waters within the Atlantic and Pacific Oceans (Chl < 0.5 mg m^{−3}). To verify the consistency of the algorithm performance, the current (SeaDAS) Iterative NIR (INIR) aerosol correction scheme [26] was applied to MODIS-Aqua data together with the refined and SeaDAS Rayleigh-scattering radiance algorithms. For typical oceanic waters, the common black-ocean assumption of negligible NIR water-leaving radiance is generally valid and proved successful for estimating aerosol radiances and generating ocean color products [25] although the INIR scheme showed significant improvements in the *Lwn* retrievals from time series SeaWiFS data.

Figure 8 shows a similar and useful comparison of *Lwn* spectra for the open oceanic waters using the refined and SeaDAS algorithms. There are significant spectral departures between the values of *Lwn* retrieved using the SeaDAS algorithm and that measured with in-situ radiometers. When retrieved using the refined algorithm, corresponding spectral plots for the entire visible wavelengths reveal no significant difference between the retrieved and measured values of *Lwn*. When compared with the NOMAD in-situ data, the accuracy of *Lwn* retrievals by the refined algorithm is improved by 12%, 13%, 15%, 6%, and 165% at 412, 443, 488, 555, and 678 nm, respectively. A significant departure of *Lwn* from the SeaDAS algorithm arises because it underestimates *L _{r}* (Δ

*L*(%) −2.358, −2.158, −1.969, −1.812, and −1.686 at 412, 443, 488, 555, and 678 nm) with an increasing degradation of the retrievals toward the shorter wavelengths and overestimates

_{r}*Lwn*(Δ

*Lwn*(%) 21.52, 19.77, 17.23, 39.56, and 79.59 at 412, 443, 488, 555, 667 and 678 nm, respectively).

Statistical comparison gives further insights into the difference between the retrieved and measured values of *Lwn* using refined and SeaDAS algorithms (Table 2). It is seen that MRE and RMSE become significantly smaller for the refined algorithm than for the SeaDAS algorithm. The correlation coefficents are also close to unity for the blue wavelengths for the refined algorithm. The slightly larger error MRE and RMSE in the retrieved *Lwn* at longer wavelengths (red and NIR) using both refined and SeaDAS algorithms is likely due to the aerosol correction. While the SeaDAS Rayleigh-scattering radiance algorithm compromises the quality of atmospheric correction and degrades the quality and consistency of derived ocean color products as observed in clear oceanic and turbid coastal waters, employing the refined algorithm yields more consistent *Lwn* retrievals (and thus more consistent biogeochemical products) and a better atmospheric connection at all visible and NIR wavelengths under different oceanic and atmospheric conditions.

#### 5.2 Algorithm evaluation using satellite data

To facilitate further analysis of the effects of *L _{r}* on satellite-derived

*Lwn*products, this section presents an application of the refined and SeaDAS algorithms to MODIS-Aqua imagery acquired over the different regional waters under different atmospheres and solar zenith angles. For this application, four MODIS-Aqua scenes were acquired on 25 October 2016 in the Pacific Ocean near Hawaii (MOBY site at around 19.5429° N, 155.6659° W), 17 March 2016 in the North Atlantic Ocean (LISCO ARONET-OC site at around 40.955° N, 73.342° W), 17 March 2011 in the Gulf of Mexico (WaveCIS AERONET-OC site at around 28.867° N, 90.483° W), and 17 January 2019 in the Southern Ocean (at around 64.09° S, 60.58° E).

Figure 9 presents examples of the MODIS-Aqua derived *L _{r}* and

*Lwn*products along with histograms as a function of the

*L*and

_{r}*Lwn*and the relative errors (Δ

*L*and Δ

_{r}*Lwn*) at the MODIS-Aqua 412 nm band for the Pacific Ocean near Hawaii (MOBY). In this case the values of

*L*at 412 nm across the scene computed using the SeaDAS algorithm are considerably lower than those computed using the refined algorithm and this deviation in

_{r}*L*resulting from the SeaDAS algorithm yielded the underestimated values of

_{r}*L*(Δ

_{r}*L*= −1.65%) and overestimated values of

_{r}*Lwn*(Δ

*Lwn*= 8.15%) at 412 nm as discussed previously. The manifestation of this effect differed in the individual bands; the lower the wavelength, the higher the deviation in

*L*; this deviation is reduced with increasing wavelength (images not shown for other MODIS-Aqua bands for brevity). The underestimated values of

_{r}*L*by the SeaDAS algorithm degraded the

_{r}*Lwn*products and significantly altered spatial structure in the backscattered radiance signal of different water components. Consistent with the previous observations, the

*L*product obtained using the refined algorithm is more reasonable for this clear oceanic site as it accounts the effect of variations in pressure, composition of air and air mass in the atmosphere and produced consistent

_{r}*Lwn*retrievals as discussed previously.

Figures 10 and 11 show similar results from the MODIS-Aqua images along with histograms as a function of the *L _{r}*,

*Lwn*and the relative errors (Δ

*L*and Δ

_{r}*Lwn*) on 17 March 2016 in the North Atlantic Ocean (LISCO ARONET-OC site) and 17 March 2011 in the Gulf of Mexico (WaveCIS AERONET-OC site). Comparing these sites, the chosen North Atlantic coastal region is more dominated by particulate and dissolved substances from the Quinnipiac, Housatonic, Norwalk and Connecticut River systems [42,43]. Similar water types exist in the northern coastal areas of the Gulf of Mexico due to the influence of multiple river systems (Mississippi, Econfina, Aucilla, St. Marks, Ochlocknee, Carrabelle, and Apalachicola), but the open-sea region of the Gulf of Mexico is generally characterized by the phytoplankton and co-related constituents [44]

For these cases, the values of *L _{r}* computed using these two algorithms are much higher and vary strongly for small to moderate values of solar zenith angles and with position across the scene (from west to east) than they were in the previous case due to the larger contribution of the molecular atmosphere and the higher range of solar zenith angles. Again, the SeaDAS algorithm underestimates

*L*and begins to overestimate

_{r}*Lwn*across these MODIS-Aqua scenes. The refined algorithm still remains consistent in computing the

*L*regardless of the atmospheric conditions and producing the expected values of

_{r}*Lwn*at 443 nm since the

*Lwn*(412 nm) signal decreases with increasing constituents’ concentrations in coastal waters of the North Atlantic (peak in histogram around 0.5 µWcm

^{−1}nm

^{−1}sr

^{−1}at 412 nm) and increases with decreasing constituents’ concentrations in relatively clear open-sea waters of the Gulf of Mexico (peak in histogram around 1 µWcm

^{−1}nm

^{−1}sr

^{−1}at 412 nm) (note the more coastal-dominated scene in the North Atlantic and the more open-sea dominated scene in the Gulf of Mexico). The relative error Δ

*L*for the wavelength 412 was −1.61% and −1.63%, which corresponded to the error Δ

_{r}*Lwn*of 33.68% and 21.15%, for the MODIS-Aqua scenes from the North Atlantic Ocean and Gulf of Mexico, respectively.

To see the effect of these computations over the southern polar ocean, the results were compared using two MODIS-Aqua scenes acquired over the Southern Ocean (17 January 2019) and Northern Ocean (not shown for brevity). In this case the solar zenith angle is larger than that observed by the MODIS-Aqua sensor over the MOBY and AERONET-OC sites. The two MODIS-Aqua scenes acquired over the polar regions were the cases with clear atmosphere and homogeneous ocean water (Fig. 12). As described earlier, similar results were obtained using the SeaDAS algorithm for the Southern Ocean with the significant departure in the computed values of *L _{r}* (Δ

*L*= −1.32 at 412 nm). The inaccurate

_{r}*L*computations by the SeaDAS algorithm caused (overestimation) errors (Δ

_{r}*Lwn*) of up to 3.70% in the retrieved values of

*Lwn*at 412 nm. These results indicate that the standard algorithm implemented in the SeaDAS system produces spatially and temporally less consistent

*L*products in comparison to the refined algorithm, which could cause significant errors in the derived

_{r}*Lwn*products at visible bands under different atmospheric and solar zenith angle conditions.

#### 5.3 The effect of the compensation factor (ROT/airmass)

The intensity of *L _{r}* is computed as a function of the ROT and airmass for various atmospheric conditions and all possible solar-sensor geometries. To take into account and accommodate the additional contribution due to variations in these parameters at different latitudes and altitudes, a compensation factor (ROT/airmass ratio) was introduced and incorporated in the computation of the

*L*using the refined algorithm. Thus, it is desirable to understand and evaluate the effect of this compensation factor on the satellite-retrieved

_{r}*L*and

_{r}*Lwn*products in various oceanic and inland water cases. For this analysis, two cases were considered: (i) a high-altitude lake Lake Titicaca which is the world’s second largest (covering 8300 square km) and highest (3810 meters above sea level) navigable lake situated in the Andes Mountains of South America between Peru and Bolivia, and (ii) many coastal and open-oceanic cases across the latitudes (MOBY, AERONET-OC and NOMAD as shown in Figs. 5 and 8). Tests carried out to determine the effect of the ROT/airmass ratio values (5%, 10%, 20% and 30% contributions considered) on the

*L*and

_{r}*Lwn*using MODIS-Aqua data showed that though quantitatively the small difference between the values of

*L*for the calculations, employing the ROT/airmass ratio above a threshold 10% begins to yield significantly low or even negative

_{r}*Lwn*values in Titicaca Lake waters with the percent error of 40-100% (Δ

*Lwn*) for blue bands (Fig. 13). Also, the magnitude of

*Lwn*retrievals using the refined algorithm with the 10% ROT/airmass effect was consistently lower than that resulting from the SeaDAS algorithm and the 5% ROT/airmass effect.

In the absence of measured in-situ data for this lake case, the MODIS-Aqua retrieved *Lwn* values are weak in the NIR wavebands and much higher in the blue wavebands for low-turbid lake waters. The spectral form and magnitude of the MODIS-Aqua derived *Lwn* are analogous to similar cases of high latitude lakes and low-turbid oceanic waters considered in this study. Tests carried out for clear (MOBY and NOMAD) and turbid (AERONET-OC) oceanic cases revealed similar results (bottom panels), where the refined algorithm with the 20% and 30% ROT/airmass effects significantly overestimated *L _{r}* and underestimated

*Lwn*for blue bands with increasing errors for turbid coastal oceanic waters (MRE 60-100%) than for open-oceanic waters (MRE up to 30%). The performance of the SeaDAS algorithm remained consistent in terms of underestimating

*L*and overestimating

_{r}*Lwn*with significant errors in the blue bands (MRE up to 70% for open oceanic cases and 280% for turbid oceanic cases). Despite the slight improvement in the

*Lwn*retrievals with the 5% ROT/airmass effect using the refined algorithm, the magnitude of the error is still significant (MRE from 50% for open-oceanic waters to 150% for turbid coastal oceanic waters). Conversely, employing the 10% ROT/airmass effect in the refined algorithm was sufficiently accurate and adequately accommodated the additional contributions due to a Rayleigh-scattering atmosphere across the latitudes and altitudes, and yielded good agreement between MODIS-Aqua and ship measurements of

*Lwn*for both oceanic and high-altitude lake waters. Significant errors in the computations of the

*L*using the SeaDAS algorithm and higher ROT/airmass effects likely affected the quality and accuracy of aerosol correction leading to much larger errors in the

_{r}*Lwn*retrievals for turbid oceanic and high-altitude lake waters. By contrast, the refined algorithm with the 10% ROT/airmass effect yielded good retrievals of

*Lwn*(and thus biogeochemical products) for diverse oceanic and atmospheric conditions. This approximation is flexible and allows for a better refinement of the algorithm based on the ROT and airmass values for accurate computation of the

*L*under different atmospheric conditions. Based on these comparisons it is believed that the refined algorithm (with the 10% RTO/airmass effect) presented here is sufficiently accurate to accommodate the variations in Rayleigh-scattering atmospheres to determine the

_{r}*L*for visible and NIR wavebands and all possible solar-sensor geometries. The computations above can easily be adapted for use in satellite data processing for oceanic and inland water applications.

_{r}## 6. Conclusion

Precise computations of the Rayleigh-scattering radiance is required to achieve desirable accuracies and spatial and temporal consistencies in the derived water-leaving radiance and biogeochemical products from satellite data because this radiance contribution to the sensor recorded TOA signal at the blue wavelengths often exceeds 80% and a few percent uncertainty in its computation can translate into several percent variation in the derived water-leaving radiance at the blue wavelengths. In this work, a refined algorithm has been developed for computations of the Rayleigh-scattering radiance and for application to the current and future water color space borne remote sensing instruments with higher radiometric sensitivity. Unlike the standard algorithm, an independent model has been developed for accurate calculations of the Rayleigh-scattering optical thickness with the effect of variations in atmospheric conditions induced by variations in surface pressure, CO_{2} volume concentration, and relative air mass for given sun-sensor geometries. Finally, the refined algorithm has been incorporated into SeaDAS system and the Rayleigh lookup tables have been regenerated with the above modifications for the computation of the spectral Rayleigh-scattering radiances for use with the water color remote sensing instruments.

For a quantitative assessment of the algorithm, several tests were carried out using MODIS-Aqua data acquired over the open ocean and coastal regions. First, the actual values of *L _{r}* and

*Lwn*were obtained using the refined algorithm for the cases where in-situ measurements are available (MOBY and AERONET-OC). Second, the actual values of

*L*and

_{r}*Lwn*were obtained for the cases where MODIS-Aqua scenes are acquired under various solar zenith angles and atmospheric conditions. The results of the

*L*computations and

_{r}*Lwn*retrievals obtained using the refined algorithm were then compared with those obtained using the SeaDAS algorithm. The results indicate that there is the strong variation of the relative percent error in the computations of the

*L*for visible bands; in particular, the error Δ

_{r}*L*induced by the SeaDAS algorithm becomes progressively larger as the wavelength gets shorter (Δ

_{r}*L*from −1.69 at 678 nm to −2.39% at 412 nm). A few percent error in the computed

_{r}*L*values eventually translates into several percent error in the retrieved

_{r}*Lwn*values (Δ

*Lwn*from 12.02% at 678 to 64.42% at 412 nm; MRE from 133.6% at 678 nm to 51.8% at 412 nm; higher MRE values at longer wavelengths due to the aerosol correction). The overall improvements achieved by the refined algorithm (relative to the SeaDAS algorithm) were 77%, 40%, 14%, 10%, and 23% for the wavelengths at 412, 443, 488, 555, and 678nm (MRE from 14% at 678 nm to 46% at 412 nm). The higher MRE values between the retrieved and measured

*Lwn*values are due to the aerosol correction. However, results indicate that significant departures between the values of

*Lwn*retrieved using the SeaDAS algorithm and that measured with field instruments are caused by the computations of

*L*values (underestimation). In contrast, the Rayleigh-scattering radiances computed using the refined algorithm produced spatially and temporally consistent water-leaving radiance retrievals as compared to the in-situ measurements regardless of the atmospheric and oceanic conditions.

_{r}Application to MODIS-Aqua images across the latitudes acquired under low-high solar zenith angles and various atmospheric conditions showed similar departures in the computed values of *L _{r}* and retrieved values of

*Lwn*. Significant errors in computations of the

*L*(underestimation) could result from the inadequacy of the SeaDAS algorithm to fully account for the variations in the atmospheric pressure, composition of air and air mass for given solar-sensor geometries. For all test cases, the refined algorithm remains stable in producing spatially and temporally consistent

_{r}*L*and

_{r}*Lwn*products for various atmospheres and solar zenith angles. With the development of new generation water color sensors having UV bands for coastal and inland applications, this problem will become more serious with the standard algorithm because of the increasing departure between the computed and actual values of

*L*at the blue and UV wavelengths. Such errors can seriously compromise the accuracy of atmospheric correction and degrade the quality of derived biogeochemical products from satellite data. Results based on the MODIS-Aqua and in-situ measurements suggest that employing the refined algorithm for computation of the

_{r}*L*can yield more accurate

_{r}*Lwn*and biogeochemical products at different latitudes and altitudes as desired by the scientific community.

## Funding

National Key Research and Development Program of China (2017YFA0603003); Natural Resources Data Management System (DST), India (OEC1819150DSTXPSHA); Indian Institute of Technology Madras.

## Acknowledgments

We would like to thank the contributors of AERONET-OC data (Bill Gibson, Robert Arnone, Sherwin Ladner, Hui Feng, Heidi M. Sosik, Dimitry Van der Zande, Giuseppe Zibordi, Sam Ahmed, Alex Gilerson) and MOBY station and NOAA Greenhouse Gas Marine Boundary Layer Reference. We are thankful to the NASA OBPG for developing and maintaining the SeaDAS software package and providing the MODIS-Aqua data and NOMAD data. We duly acknowledge the MOBY operations team (Hawaii ﬂoats) for generating and making the radiometric data available for research purpose. We sincerely thank the two anonymous reviewers for their valuable comments and suggestions which helped to improve the quality of this manuscript.

## References

**1. **R. K. Singh and P. Shanmugam, “Corrigendum to “A novel method for estimation of aerosol radiance and its extrapolation in the atmospheric correction of satellite data over optically complex oceanic waters” [Remote Sensing of Environment 142 (2014) 188–206],” Remote Sens. Environ. **148**, 222–223 (2014). [CrossRef]

**2. **R. K. Singh and P. Shanmugam, “A novel method for estimation of aerosol radiance and its extrapolation in the atmospheric correction of satellite data over optically complex oceanic waters,” Remote Sens. Environ. **142**, 188–206 (2014). [CrossRef]

**3. **P. Shanmugam, “A new bio-optical algorithm for the remote sensing of algal blooms in complex ocean waters,” J. Geophys. Res. **116**(C4), C04016 (2011). [CrossRef]

**4. **C. Schueler, J. Yoder, D. Antoine, C. Castillo, R. Evans, C. Mengelt, C. Mobley, J. Sarmiento, S. Sathyendranath, D. Siegel, and C. Wilson, “Assessing Requirements for Sustained Ocean Color Research and Observations,” in AIAA SPACE 2011 Conference & Exposition (American Institute of Aeronautics and Astronautics, 2011), pp. 1–7.

**5. **K. Ruddick, Q. Vanhellemont, J. Yan, G. Neukermans, G. Wei, and S. Shang, “Variability of suspended particulate matter in the Bohai Sea from the geostationary Ocean Color Imager (GOCI),” Ocean Sci. J. **47**(3), 331–345 (2012). [CrossRef]

**6. **J. Zhao, M. Temimi, and H. Ghedira, “Characterization of harmful algal blooms (HABs) in the Arabian Gulf and the Sea of Oman using MERIS fluorescence data,” ISPRS J. Photogramm. Remote Sens. **101**, 125–136 (2015). [CrossRef]

**7. **S. Sathyendranath, F. E. Hoge, T. Platt, and R. N. Swift, “Detection of phytoplankton pigments from ocean color: improved algorithms,” Appl. Opt. **33**(6), 1081–1089 (1994). [CrossRef]

**8. **N. Pahlevan, J. R. Schott, B. A. Franz, G. Zibordi, B. Markham, S. Bailey, C. B. Schaaf, M. Ondrusek, S. Greb, and C. M. Strait, “Landsat 8 remote sensing reflectance (Rrs) products: Evaluations, intercomparisons, and enhancements,” Remote Sens. Environ. **190**, 289–301 (2017). [CrossRef]

**9. **J. E. O’Reilly, S. Maritorena, B. G. Mitchell, D. A. Siegel, K. L. Carder, S. A. Garver, M. Kahru, and C. McClain, “Ocean color chlorophyll algorithms for SeaWiFS,” J. Geophys. Res. **103**(C11), 24937–24953 (1998). [CrossRef]

**10. **F. Gerald Kelly, in * Minimum Requirements for an Operational, Ocean-Color Sensor for the Open Ocean* (International Ocean-Color Coordinating Group, 1998).

**11. **C. Hu, K. L. Carder, and F. E. Muller-Karger, “Atmospheric Correction of SeaWiFS Imagery over Turbid Coastal Waters,” Remote Sens. Environ. **74**(2), 195–206 (2000). [CrossRef]

**12. **Y.-H. Ahn and P. Shanmugam, “Detecting the red tide algal blooms from satellite ocean color observations in optically complex Northeast-Asia Coastal waters,” Remote Sens. Environ. **103**(4), 419–437 (2006). [CrossRef]

**13. **S. W. Bailey and P. J. Werdell, “A multi-sensor approach for the on-orbit validation of ocean color satellite data products,” Remote Sens. Environ. **102**(1-2), 12–23 (2006). [CrossRef]

**14. **P. Shanmugam, “New models for retrieving and partitioning the colored dissolved organic matter in the global ocean: Implications for remote sensing,” Remote Sens. Environ. **115**(6), 1501–1521 (2011). [CrossRef]

**15. **H. R. Gordon, “Atmospheric correction of ocean color imagery in the Earth Observing System era,” J. Geophys. Res.: Atmos. **102**(D14), 17081–17106 (1997). [CrossRef]

**16. **H. R. Gordon, “In-Orbit Calibration Strategy for Ocean Color Sensors,” Remote Sens. Environ. **63**(3), 265–278 (1998). [CrossRef]

**17. **M. Wang, “A refinement for the Rayleigh radiance computation with variation of the atmospheric pressure,” Int. J. Remote Sens. **26**(24), 5651–5663 (2005). [CrossRef]

**18. **C. Tomasi, V. Vitale, B. Petkov, A. Lupi, and A. Cacciari, “Improved algorithm for calculations of Rayleigh-scattering optical depth in standard atmospheres,” Appl. Opt. **44**(16), 3320–3341 (2005). [CrossRef]

**19. **H. R. Gordon, J. W. Brown, and R. H. Evans, “Exact Rayleigh scattering calculations for use with the Nimbus-7 Coastal Zone Color Scanner,” Appl. Opt. **27**(5), 862–871 (1988). [CrossRef]

**20. **H. R. Gordon and M. Wang, “Surface-roughness considerations for atmospheric correction of ocean color sensors II: Error in the retrieved water-leaving radiance,” Appl. Opt. **31**(21), 4247–4260 (1992). [CrossRef]

**21. **P. M. Teillet, “Rayleigh optical depth comparisons from various sources,” Appl. Opt. **29**(13), 1897–1900 (1990). [CrossRef]

**22. **L. Vinet and A. Zhedanov, “A 'missing' family of classical orthogonal polynomials,” Space Sci. Rev. **16**, 527–610 (2010). [CrossRef]

**23. **A. J. Brown, “Spectral bluing induced by small particles under the Mie and Rayleigh regimes,” Icarus **239**, 85–95 (2014). [CrossRef]

**24. **G. Zibordi, F. Mélin, J.-F. Berthon, B. Holben, I. Slutsker, D. Giles, D. D’Alimonte, D. Vandemark, H. Feng, G. Schuster, B. E. Fabbri, S. Kaitala, and J. Seppälä, “AERONET-OC: A Network for the Validation of Ocean Color Primary Products,” J. Atmos. Ocean. Technol. **26**(8), 1634–1651 (2009). [CrossRef]

**25. **H. R. Gordon and M. Wang, “Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm,” Appl. Opt. **33**(3), 443–452 (1994). [CrossRef]

**26. **S. W. Bailey, B. A. Franz, and P. J. Werdell, “Estimation of near-infrared water-leaving reflectance for satellite ocean color data processing,” Opt. Express **18**(7), 7521–7527 (2010). [CrossRef]

**27. **M. Wang, S. Son, and W. Shi, “Evaluation of MODIS SWIR and NIR-SWIR atmospheric correction algorithms using SeaBASS data,” Remote Sens. Environ. **113**(3), 635–644 (2009). [CrossRef]

**28. **T. Varunan and P. Shanmugam, “An optical tool for quantitative assessment of phycocyanin pigment concentration in cyanobacterial blooms within inland and marine environments,” J. Great Lakes Res. **43**(1), 32–49 (2017). [CrossRef]

**29. **H. J. Nasiha, P. Shanmugam, and R. Sundaravadivelu, “Estimation of sediment settling velocity in estuarine and coastal waters using optical remote sensing data,” Adv. Space Res. **63**(11), 3473–3488 (2019). [CrossRef]

**30. **T. Varunan and P. Shanmugam, “Use of Landsat 8 data for characterizing dynamic changes in physical and acoustical properties of coastal lagoon and estuarine waters,” Adv. Space Res. **62**(9), 2393–2417 (2018). [CrossRef]

**31. **B. A. Bodhaine, N. B. Wood, E. G. Dutton, and J. R. Slusser, “On Rayleigh Optical Depth Calculations,” J. Atmos. Ocean. Technol. **16**(11), 1854–1861 (1999). [CrossRef]

**32. **R. Penndorf, “Tables of the Refractive Index for Standard Air and the Rayleigh Scattering Coefficient for the Spectral Region between 02 and 200 μ and Their Application to Atmospheric Optics,” J. Opt. Soc. Am. **47**(2), 176–182 (1957). [CrossRef]

**33. **E. G. Dutton, P. Reddy, S. Ryan, and J. J. DeLuisi, “Features and effects of aerosol optical depth observed at Mauna Loa, Hawaii: 1982–1992,” J. Geophys. Res. **99**(D4), 8295–8306 (1994). [CrossRef]

**34. **M. Nicolet, “On the molecular scattering in the terrestrial atmosphere : An empirical formula for its calculation in the homosphere,” Planet. Space Sci. **32**(11), 1467–1468 (1984). [CrossRef]

**35. **H. R. Gordon and M. Wang, “Surface-roughness considerations for atmospheric correction of ocean color sensors 1: The Rayleigh-scattering component,” Appl. Opt. **31**(21), 4247–4260 (1992). [CrossRef]

**36. **B. Edlén, “The Refractive Index of Air,” Metrologia **2**(2), 71–80 (1966). [CrossRef]

**37. **D. R. Bates, “Rayleigh scattering by air,” Planet. Space Sci. **32**(6), 785–790 (1984). [CrossRef]

**38. **R. J. List, “Smithsonian meteorological tables,” Smithsonian Miscellaneous Collections **114**(1), 1–527 (1968).

**39. **C. Tomasi and B. H. Petkov, “Spectral calculations of Rayleigh-scattering optical depth at Arctic and Antarctic sites using a two-term algorithm,” J. Geophys. Res.: Atmos. **120**(18), 9514–9538 (2015). [CrossRef]

**40. **L. W. Thomason, B. M. Herman, and J. A. Reagan, “The Effect of Atmospheric Attenuators with Structured Vertical Distributions on Air Mass Determinations and Langley Plot Analyses,” J. Atmos. Sci. **40**(7), 1851–1854 (1983). [CrossRef]

**41. **F. Kasten and A. T. Young, “Revised optical air mass tables and approximation formula,” Appl. Opt. **28**(22), 4735–4738 (1989). [CrossRef]

**42. **M. Meybeck, “Atmospheric inputs and river transport of dissolved substances. In: Dissolved Loads of 0.Rivers and Surface Water Quantity/Quality Relationships,” in Proceedings of the Hamburg Symposium (IAHS Publication, 1983), pp. 173–192.

**43. **P. A. Raymond and J. E. Bauer, “Riverine export of aged terrestrial organic matter to the North Atlantic Ocean,” Nature **409**(6819), 497–500 (2001). [CrossRef]

**44. **F. A. C. Jonathan R, Joseph N. Pennock, Jorge A. Boyer, Richard L. Herrern-Silveira, Terry E. Iverson, and Behzad Mortazavi Whitledgc, “Nutrient Behavior and Phytoplankton Production in Gulf of Mexico Estuaries,” in * Biogeochemistry of Gulf of Mexico Estuaries*, R. R. T. Thomas, S. Bianchi, and Jonathan R. Pennock, eds. (John Wiley & Sons, Inc., 1999), pp. 109–162.