Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

New algorithm for computation of the Rayleigh-scattering radiance for remote sensing of water color from space

Vimalathitthan Shanmugam, Palanisamy Shanmugam, and Xianqiang He

Open Access Open Access

Abstract

Rayleigh-scattering radiance (Lr) 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 (CO2 volume concentration), surface atmospheric pressure and relative air mass for given sun-sensor geometries. As CO2 is well mixed throughout the atmospheric column, the CO2 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 Lr, particularly in the UV-Blue region of the spectrum. The refined algorithm for calculating the Lr is tested on several MODIS-Aqua Level 1A data and the relative errors in Rayleigh-scattering radiance and normalized water-leaving radiance (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 Lr calculated using the SeaDAS algorithm are mostly underestimated and show significant departures with the Lr calculated using the refined algorithm. This departure induced by the SeaDAS algorithm to Lr becomes larger with decreasing wavelength (ΔLr from −2.38% at 412 nm to 1.69% at 678 nm), which causes errors in 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 Lr can yield more accurate 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, Lr), 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 O2, N2, CO2, 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 Lr 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 Lr 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 Lr 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 Lr calculations for visible bands.

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 N2, O2, Ar, CO2 and H2O) [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 Lr 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, CO2 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 Lr and 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.

 figure: Fig. 1.

Fig. 1. Map of in-situ sampling location points in coastal and open ocean waters.

Download Full Size | PDF

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_{TOA}}(\lambda ) = {L_r}(\lambda ) + {L_a}(\lambda ) + {L_{ra}}(\lambda ) + T(\lambda ){L_g}(\lambda ) + {L_{wc}}(\lambda ) + t(\lambda ){L_w}(\lambda ).$$
where LTOA(λ), Lr(λ), La(λ), Lra(λ), Lg(λ), Lwc(λ) and Lw(λ) 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]. 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 Lg(λ) and Lwc(λ) under high wind speed conditions, these contributions are generally ignored for brevity. The contributions La(λ) and Lra(λ) 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].

Because Lr(λ) is the most dominant signal in the LTOA(λ), it is critical that the Rayleigh-scattering radiance must be computed accurately. According to Gordon et al. [19], Lr(λ) was computed using the single scattering approximation for analysis of CZCS imagery,

$${L_r}(\lambda ) = {{{F^{\prime}}_0(\lambda ){\tau _r}(\lambda ){p_r}({\theta _0},\theta ,\lambda )} \mathord{\left/ {\vphantom {{{F^{\prime}}_0(\lambda ){\tau_r}(\lambda ){p_r}({\theta_0},\theta ,\lambda )} {4\pi }}} \right.} {4\pi }}.$$
where
$${p_r}({{\theta_0},\theta ,\lambda } )= {{\{{{{\mathop{\rm P}\nolimits}_r}({{\theta_{0 - }},\lambda } )+ [{\rho ({\theta_0}) + \rho (\theta )} ]{{\mathop{\rm P}\nolimits}_r}({{\theta_{0 + }},\lambda } )} \}} \mathord{\left/ {\vphantom {{\{{{{\mathop{\rm P}\nolimits}_r}({{\theta_{0 - }},\lambda } )+ [{\rho ({\theta_0}) + \rho (\theta )} ]{{\mathop{\rm P}\nolimits}_r}({{\theta_{0 + }},\lambda } )} \}} {\cos {\theta_0}}}} \right.} {\cos {\theta _0}}}.$$
$$\cos {\theta _ \pm } = \pm \cos \theta \cos {\theta _0} - \sin \theta \sin {\theta _0}\cos ({\vartheta _0} - \vartheta ).$$
where F’0 = instantaneous extraterrestrial solar irradiance, τr = Rayleigh-scattering optical thickness, ρ(θ0) = Fresnel reflectance of the interface for an incident angle θ0, and Pr (α, λ) = Rayleigh scattering phase function. θ, ϑ, θ0, and ϑ0, are the zenith and azimuth angles of a vector from the sea surface to the sun and sensor, respectively [19]. Pr(α) is given by,
$${{\mathop{\rm P}\nolimits} _r}(\alpha ) = {3 \mathord{\left/ {\vphantom {3 4}} \right.} 4}[1 + {\cos ^2}\alpha ].$$
In Eq. (2), the most important quantity is the optical thickness of the atmosphere which is derived from the Beer’s Law
$$\frac{{I(\lambda )}}{{{I_0}(\lambda )}} = {e^{\left( { - \frac{{\tau (\lambda )}}{{\cos \theta }}} \right)}}.$$
where I0(λ) is the extraterrestrial solar radiation flux at wavelength λ, 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
$${\tau _0}(\lambda ) = 0.008569{\lambda ^{ - 4}}(1 + 0.0113{\lambda ^{ - 2}} + 0.00013{\lambda ^{ - 4}}).$$
where τo (λ) is the ROT using the standard atmosphere pressure. For the atmospheric pressure variation at a site, ROT is defined as
$${\tau _R}(\lambda ) = {\tau _0}(\lambda )\frac{P}{{{P_0}}}.$$
where P = site-specific atmospheric pressure, and P0 = 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
$$\frac{{{L_r}[{{\tau_r}({\lambda ,{P_0} + \nabla P} )} ]}}{{{L_r}[{{\tau_r}({\lambda ,{P_0}} )} ]}} = \frac{{1 - \exp [{ - {\tau_r}({\lambda ,{P_0} + \nabla P} )/\cos \theta } ]}}{{1 - \exp [{ - {\tau_r}({\lambda ,{P_0}} )/\cos \theta } ]}}.$$
where θ is the zenith angle, τr (λ, P0+∇P) and τr (λ, P) are the Rayleigh optical thicknesses, at a given wavelength, for atmospheric pressures of P0+∇P and P0 respectively. Later, Wang [17] redefined the Gordon [19] algorithm as
$$\frac{{{L_r}[{{\tau_r}({\lambda ,{P_0} + \nabla P} )} ]}}{{{L_r}[{{\tau_r}({\lambda ,{P_0}} )} ]}} = \frac{{1 - \exp [{ - C(\lambda ,M){\tau_r}({\lambda ,{P_0} + \nabla P} )M} ]}}{{1 - \exp [{ - C(\lambda ,M){\tau_r}({\lambda ,{P_0}} )M} ]}}.$$
where M = 1/cos θ0+1/cos θ is the air mass value for the solar and sensor zenith angles of θ0 and θ, 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 (Lr) is computed for a given wavelength by accounting for the variations of the CO2 volume concentration, pressure and air mass within the atmospheric column. Figure 2 describes the implementation of this scheme for Lr computations. The refined formula takes the form

$${L_r}({\tau _r},P) = {L_r}({{\tau_0}(\lambda ,{P_0})} )\left( {\frac{P}{{{P_0}}} + C} \right).$$
where P is the pressure at a site, P0 is the standard pressure (1013.25 mb), and Lr (τ0(λ, P0)) is the Rayleigh-scattering radiance for the standard atmosphere. C is the dynamic constant derived from the air mass and optical thickness of the atmosphere as
$$C = 0.1 \times \left( {\frac{{{\tau_r}}}{M}} \right).$$
where τr is the ROT calculated as a function of wavelength and CO2 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 τr 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
$$\sigma = \frac{{24{\pi ^3}{{(n_s^2 - 1)}^2}}}{{{\lambda ^4}N_s^2{{(n_s^2 + 2)}^2}}}\left( {\frac{{6 + 3\rho }}{{6 - 7\rho }}} \right).$$
where σ is the scattering cross section per molecule (cm2 molecule−1), Ns is molecular number density per unit volume, 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, Ns depends on temperature and pressure, and the resulting σ depends on the composition of the gas but is independent of temperature and pressure [31].

 figure: Fig. 2.

Fig. 2. Flowchart depicting the implementation of refined algorithm for computing the Rayleigh-scattering radiance.

Download Full Size | PDF

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

$$({n_{300}} - 1) = 8060.51 + \frac{{2480990}}{{132.274 - {\lambda ^2}}} + \frac{{17455.7}}{{39.32957 - {\lambda ^2}}}.$$
The variability in refractive index due to CO2 is given by
$$\frac{{({n_{C{O_2}}} - 1)}}{{({n_{300}} - 1)}} = 1 + 0.54(C{O_2} - 0.0003).$$
The depolarization of air is estimated as a function of the atmospheric constituents such as O2, N2, Ar and CO2 [37] as
$$F(air,C{O_2}) = \frac{{(78.084 \times x \times F({N_2}) + 20.946 \times x \times F({O_2}) + 0.934 \times x \times {C_{C{O_2}}} \times 1.15)}}{{100}}.$$
where the concentration of CO2 (CCO2) is expressed in parts per volume by percent, and x is a correction factor for other gas concentrations (according to CO2 concentration). The depolarization of N2, O2, Ar and CO2 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
$$\tau (\lambda ,C{O_2}) = \sigma \frac{{PA}}{{{m_a}g}}.$$
where P = pressure, A = Avogadro’s number, ma = mean molecular weight of the air, and 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 CO2 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 CO2 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 CO2 concentration. Latitude-wise monthly mean atmospheric CO2 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 CO2 concentration. The monthly average increment of CO2 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 CO2 concentrations and the corresponding atmospheric CO2 concentrations calculated from the fitting equations.

 figure: Fig. 3.

Fig. 3. Comparison of the monthly NOAA mean CO2 concentration and modelled CO2 concentration as a function of the latitude for the periods 2003 and 2004.

Download Full Size | PDF

Tables Icon

Table 1. The average increment per year of CO2 for each month calculated from 2000 to 2006 (NOAA Greenhouse Gas Marine Boundary Layer Reference).

View Table | View all tables in this article

The results of the calculations clearly showed the strong seasonal variation in CO2 with latitude and the fitting equation estimated CO2 are in good agreement with monthly mean CO2 measurements. Based on the CO2 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,

$$M(\theta ) = \frac{1}{{{P_0}}}\sum\limits_{{Z_0}}^{{Z_\alpha }} {\frac{{\nabla P}}{{\sqrt {1 - {{\left( {\frac{{{n_0}}}{{{n_z}}}} \right)}^2}{{\left( {\frac{{{z_0}}}{z}} \right)}^2}{{\sin }^2}\theta } }}} .$$
where P0 is the standard atmospheric pressure, and n0 and nz are the refractive indices of the lower and upper atmospheric layers, respectively. z0 and 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 Lr are performed as accurately as possible because a few percent error in Lr can translate into several percent error in the retrieved water-leaving radiance products.

 figure: Fig. 4.

Fig. 4. The variation of air mass as a function of the solar and sensor zenith angles.

Download Full Size | PDF

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 F0(λ)– Extra-terrestrial solar irradiance; Lw(λ)– Water leaving radiance; Eds(λ) – Downwelling solar irradiance.

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 CO2 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 Lr 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 Lr 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 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 Lr 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 Lwn values at the short wavelengths. In contrast, the more accurate computations of Lr by the refined algorithm significantly reduced the discrepancy and produced good agreement between MODIS-Aqua and in-situ measurements of Lwn under various atmospheric and oceanic conditions.

 figure: Fig. 5.

Fig. 5. Comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) spectra retrieved using the refined and SeaDAS algorithms with in-situ measurements in open ocean (MOBY) and coastal water sites (AERONET-OC sites – Thornton_C-power, WaveCIS, and Helsinki). Results from the standard schemes for calculating Rayleigh-scattering radiance and aerosol radiance (SeaDAS-INIR) are also included for clear oceanic waters (MOBY) only for comparison purposes.

Download Full Size | PDF

Statistical evaluation of the data based on the relative errors in Rayleigh-scattering radiance (ΔLr = LrSeaDAS- LrRefined/LrRefined) and normalized water-leaving radiance (ΔLwn = LwnSeaDAS- LwnRefined/LwnRefined) and the mean relative error in normalized water-leaving radiance (MRE = LwnSatellite- Lwnin-situ/Lwnin-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 Lr for visible bands between the SeaDAS and refined algorithms, the error ΔLr 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 Lr eventually translated into several percent error in the derived Lwn. Considering these twelve samples (Fig. 5), the error ΔLr 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 Δ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 Lr (SeaDAS algorithm) seriously compromise the accuracy of atmospheric correction and degrade the quality of derived biogeochemical products from satellite data.

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 ΔLr 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 Lr by the SeaDAS algorithm. In this case, the error ΔLwn, corresponding to the error ΔLr, 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 Lr values (Fig. 6).

 figure: Fig. 6.

Fig. 6. Comparison of the mean relative error (MRE) in the MODIS-Aqua derived normalized water-leaving radiance (Lwn) using the SeaDAS and refined algorithms for some key wavelengths (412, 443, 488, 555, and 678 nm) at the selected MOBY and AERONET sites.

Download Full Size | PDF

The significant error in the computations of the Lr 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.

Tables Icon

Table 2. Statistical comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) using the refined and SeaDAS algorithms and the MOBY/AERONET-OC in-situ data from mostly coastal oceanic water and NOMAD in-situ data from open oceanic waters. Results are shown only for some key wavelengths and errors indicate both standard and random errors along with the correlation coefficients derived from the analysis between the retrieved and measured in-situ Lwn data.

View Table | View all tables in this article

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.

 figure: Fig. 7.

Fig. 7. Comparison of histograms in mean relative error (MRE) for retrieving the normalized water-leaving radiance (Lwn) from MODIS-Aqua images using the refined and SeaDAS algorithms for some key wavelengths (412, 443, 488, 555, and 678 nm) at the sleected MOBY and AERONET sites.

Download Full Size | PDF

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 Lr (ΔLr (%) −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 Lwn (ΔLwn (%) 21.52, 19.77, 17.23, 39.56, and 79.59 at 412, 443, 488, 555, 667 and 678 nm, respectively).

 figure: Fig. 8.

Fig. 8. Comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) spectra retrieved using the refined and SeaDAS algorithms with in-situ measurements (NOMAD) in open ocean waters within Atlantic and Paficic Oceans.

Download Full Size | PDF

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 Lr 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 Lr and Lwn products along with histograms as a function of the Lr and Lwn and the relative errors (ΔLr and ΔLwn) at the MODIS-Aqua 412 nm band for the Pacific Ocean near Hawaii (MOBY). In this case the values of Lr 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 Lr resulting from the SeaDAS algorithm yielded the underestimated values of Lr (ΔLr = −1.65%) and overestimated values of 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 Lr; this deviation is reduced with increasing wavelength (images not shown for other MODIS-Aqua bands for brevity). The underestimated values of Lr by the SeaDAS algorithm degraded the Lwn products and significantly altered spatial structure in the backscattered radiance signal of different water components. Consistent with the previous observations, the Lr 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 Lwn retrievals as discussed previously.

 figure: Fig. 9.

Fig. 9. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (25 October 2016) using the SeaDAS and refined algorithms in clear oceanic water at the MOBY site (Pacific Ocean near Hawaii) Figs. 9(a)–9(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

Download Full Size | PDF

Figures 10 and 11 show similar results from the MODIS-Aqua images along with histograms as a function of the Lr, Lwn and the relative errors (ΔLr and Δ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]

 figure: Fig. 10.

Fig. 10. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 March 2016) using the SeaDAS and refined algorithms in coastal waters in the North Atlantic Ocean (LISCO ARONET-OC site) Figs. 10(a)–10(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 March 2011) using the SeaDAS and refined algorithms in coastal and open-sea waters in the Gulf of Mexico (WaveCIS AERONET-OC site) Figs. 11(a)–11(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

Download Full Size | PDF

For these cases, the values of Lr 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 Lr and begins to overestimate Lwn across these MODIS-Aqua scenes. The refined algorithm still remains consistent in computing the Lr regardless of the atmospheric conditions and producing the expected values of 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−1nm−1sr−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−1nm−1sr−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 ΔLr for the wavelength 412 was −1.61% and −1.63%, which corresponded to the error Δ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 Lr (ΔLr = −1.32 at 412 nm). The inaccurate Lr computations by the SeaDAS algorithm caused (overestimation) errors (Δ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 Lr products in comparison to the refined algorithm, which could cause significant errors in the derived Lwn products at visible bands under different atmospheric and solar zenith angle conditions.

 figure: Fig. 12.

Fig. 12. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 January 2019) using the SeaDAS and refined algorithms in the Southern Ocean) Figs. 12(a)–12(d). For comparison purpose, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

Download Full Size | PDF

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

The intensity of Lr 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 Lr using the refined algorithm. Thus, it is desirable to understand and evaluate the effect of this compensation factor on the satellite-retrieved Lr and 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 Lr and Lwn using MODIS-Aqua data showed that though quantitatively the small difference between the values of Lr for the calculations, employing the ROT/airmass ratio above a threshold 10% begins to yield significantly low or even negative 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.

 figure: Fig. 13.

Fig. 13. Results for evaluating the effect of ROT/airmass ratios (5%, 10%, 20% And 30%) with the refined algorithm on the Lr and Lwn products in comparison with the results derived from the SeaDAS algorithm using MODIS-Aqua data from Titicaca lake waters and MODIS-Aqua and in-situ matchup data from oceanic cases (MOBY, NOMAD and AERONET-OC).

Download Full Size | PDF

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 Lr 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 Lr and overestimating 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 Lr 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 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 Lr 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 Lr 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.

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, CO2 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 Lr 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 Lr and Lwn were obtained for the cases where MODIS-Aqua scenes are acquired under various solar zenith angles and atmospheric conditions. The results of the Lr computations and 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 Lr for visible bands; in particular, the error ΔLr induced by the SeaDAS algorithm becomes progressively larger as the wavelength gets shorter (ΔLr from −1.69 at 678 nm to −2.39% at 412 nm). A few percent error in the computed Lr values eventually translates into several percent error in the retrieved 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 Lr 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.

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 Lr and retrieved values of Lwn. Significant errors in computations of the Lr (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 Lr and 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 Lr 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 Lr can yield more accurate 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 floats) 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.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)
Fig. 1.
Fig. 1. Map of in-situ sampling location points in coastal and open ocean waters.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Flowchart depicting the implementation of refined algorithm for computing the Rayleigh-scattering radiance.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Comparison of the monthly NOAA mean CO2 concentration and modelled CO2 concentration as a function of the latitude for the periods 2003 and 2004.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. The variation of air mass as a function of the solar and sensor zenith angles.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. Comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) spectra retrieved using the refined and SeaDAS algorithms with in-situ measurements in open ocean (MOBY) and coastal water sites (AERONET-OC sites – Thornton_C-power, WaveCIS, and Helsinki). Results from the standard schemes for calculating Rayleigh-scattering radiance and aerosol radiance (SeaDAS-INIR) are also included for clear oceanic waters (MOBY) only for comparison purposes.

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Comparison of the mean relative error (MRE) in the MODIS-Aqua derived normalized water-leaving radiance (Lwn) using the SeaDAS and refined algorithms for some key wavelengths (412, 443, 488, 555, and 678 nm) at the selected MOBY and AERONET sites.

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Comparison of histograms in mean relative error (MRE) for retrieving the normalized water-leaving radiance (Lwn) from MODIS-Aqua images using the refined and SeaDAS algorithms for some key wavelengths (412, 443, 488, 555, and 678 nm) at the sleected MOBY and AERONET sites.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. Comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) spectra retrieved using the refined and SeaDAS algorithms with in-situ measurements (NOMAD) in open ocean waters within Atlantic and Paficic Oceans.

View in Article | Download Full Size | PDF

Fig. 9.
Fig. 9. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (25 October 2016) using the SeaDAS and refined algorithms in clear oceanic water at the MOBY site (Pacific Ocean near Hawaii) Figs. 9(a)–9(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

View in Article | Download Full Size | PDF

Fig. 10.
Fig. 10. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 March 2016) using the SeaDAS and refined algorithms in coastal waters in the North Atlantic Ocean (LISCO ARONET-OC site) Figs. 10(a)–10(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

View in Article | Download Full Size | PDF

Fig. 11.
Fig. 11. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 March 2011) using the SeaDAS and refined algorithms in coastal and open-sea waters in the Gulf of Mexico (WaveCIS AERONET-OC site) Figs. 11(a)–11(d). For comparison purposes, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

View in Article | Download Full Size | PDF

Fig. 12.
Fig. 12. Comparison of the spatial products of Rayleigh-scattering radiance (Lr) and normalized water-leaving radiance (Lwn) derived from MODIS-Aqua data (17 January 2019) using the SeaDAS and refined algorithms in the Southern Ocean) Figs. 12(a)–12(d). For comparison purpose, histogram representations of Lr, Lwn, ΔLr and ΔLwn products derived using the SeaDAS and refined algorithms for the entire MODIS-Aqua scene are shown in the bottom panels.

View in Article | Download Full Size | PDF

Fig. 13.
Fig. 13. Results for evaluating the effect of ROT/airmass ratios (5%, 10%, 20% And 30%) with the refined algorithm on the Lr and Lwn products in comparison with the results derived from the SeaDAS algorithm using MODIS-Aqua data from Titicaca lake waters and MODIS-Aqua and in-situ matchup data from oceanic cases (MOBY, NOMAD and AERONET-OC).

View in Article | Download Full Size | PDF

Tables (2)

Tables Icon

Table 1. The average increment per year of CO2 for each month calculated from 2000 to 2006 (NOAA Greenhouse Gas Marine Boundary Layer Reference).

View Table | View all tables in this article
Tables Icon

Table 2. Statistical comparison of the MODIS-Aqua derived normalized water-leaving radiance (Lwn) using the refined and SeaDAS algorithms and the MOBY/AERONET-OC in-situ data from mostly coastal oceanic water and NOMAD in-situ data from open oceanic waters. Results are shown only for some key wavelengths and errors indicate both standard and random errors along with the correlation coefficients derived from the analysis between the retrieved and measured in-situ Lwn data.

View Table | View all tables in this article

Equations (18)

Equations on this page are rendered with MathJax. Learn more.

L T O A ( λ ) = L r ( λ ) + L a ( λ ) + L r a ( λ ) + T ( λ ) L g ( λ ) + L w c ( λ ) + t ( λ ) L w ( λ ) .
L r ( λ ) = F ′ 0 ( λ ) τ r ( λ ) p r ( θ 0 , θ , λ ) / F ′ 0 ( λ ) τ r ( λ ) p r ( θ 0 , θ , λ ) 4 π 4 π .
p r ( θ 0 , θ , λ ) = { P r ( θ 0 − , λ ) + [ ρ ( θ 0 ) + ρ ( θ ) ] P r ( θ 0 + , λ ) } / { P r ( θ 0 − , λ ) + [ ρ ( θ 0 ) + ρ ( θ ) ] P r ( θ 0 + , λ ) } cos ⁡ θ 0 cos ⁡ θ 0 .
cos ⁡ θ ± = ± cos ⁡ θ cos ⁡ θ 0 − sin ⁡ θ sin ⁡ θ 0 cos ⁡ ( ϑ 0 − ϑ ) .
P r ( α ) = 3 / 3 4 4 [ 1 + cos 2 α ] .
I ( λ ) I 0 ( λ ) = e ( − τ ( λ ) cos ⁡ θ ) .
τ 0 ( λ ) = 0.008569 λ − 4 ( 1 + 0.0113 λ − 2 + 0.00013 λ − 4 ) .
τ R ( λ ) = τ 0 ( λ ) P P 0 .
L r [ τ r ( λ , P 0 + ∇ P ) ] L r [ τ r ( λ , P 0 ) ] = 1 − exp ⁡ [ − τ r ( λ , P 0 + ∇ P ) / cos ⁡ θ ] 1 − exp ⁡ [ − τ r ( λ , P 0 ) / cos ⁡ θ ] .
L r [ τ r ( λ , P 0 + ∇ P ) ] L r [ τ r ( λ , P 0 ) ] = 1 − exp ⁡ [ − C ( λ , M ) τ r ( λ , P 0 + ∇ P ) M ] 1 − exp ⁡ [ − C ( λ , M ) τ r ( λ , P 0 ) M ] .
L r ( τ r , P ) = L r ( τ 0 ( λ , P 0 ) ) ( P P 0 + C ) .
C = 0.1 × ( τ r M ) .
σ = 24 π 3 ( n s 2 − 1 ) 2 λ 4 N s 2 ( n s 2 + 2 ) 2 ( 6 + 3 ρ 6 − 7 ρ ) .
( n 300 − 1 ) = 8060.51 + 2480990 132.274 − λ 2 + 17455.7 39.32957 − λ 2 .
( n C O 2 − 1 ) ( n 300 − 1 ) = 1 + 0.54 ( C O 2 − 0.0003 ) .
F ( a i r , C O 2 ) = ( 78.084 × x × F ( N 2 ) + 20.946 × x × F ( O 2 ) + 0.934 × x × C C O 2 × 1.15 ) 100 .
τ ( λ , C O 2 ) = σ P A m a g .
M ( θ ) = 1 P 0 ∑ Z 0 Z α ∇ P 1 − ( n 0 n z ) 2 ( z 0 z ) 2 sin 2 θ .
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use