## Abstract

The atmospheric correction algorithm employed by the NASA Ocean Biology Processing Group requires an assumption of negligible water-leaving reflectance in the near-infrared region of the spectrum. For waters where this assumption is not valid, an optical model is used to estimate near-infrared water-leaving reflectance. We describe this optical model as implemented for the sixth reprocessing of the SeaWiFS mission-long time-series (September 2009). Application of the optical model resulted in significant reductions in the number of negative water-leaving reflectance retrievals in turbid and optically complex waters, and improved agreement with in situ chlorophyll-a observations. The incidence of negative water-leaving reflectance retrievals at 412 nm was reduced by 40%, while negative reflectance at 490 nm was nearly eliminated.

©2010 Optical Society of America

## 1. Introduction

The retrieval of accurate geophysical data products (e.g., spectral remote sensing reflectance R_{rs}(λ), and concentrations of the phytoplankton pigment chlorophyll-a, C_{a}) from satellite-borne radiometers requires the use of an ‘atmospheric correction’ algorithm [1]. A critical component of this atmospheric correction algorithm is estimation of spectral aerosol reflectance, ρ_{a}(λ), which are subtracted from the total measured signal as part of the process of determining R_{rs}(λ). The algorithm implemented by the NASA Ocean Biology Processing Group (OBPG) [1] requires an assumption of negligible water-leaving reflectance in the near infrared (NIR) region of the spectrum (i.e., R_{rs}(NIR) = 0 sr^{−1}). The NIR bands of the NASA Sea-viewing Wide Field-of-view Sensor (SeaWiFS) are centered on 765 and 865 nm. With this ‘black-pixel’ assumption, the measured top-of-atmosphere (TOA) reflectance in two NIR bands can be used to estimate both the magnitude and spectral dependence of ρ_{a}(NIR). Unfortunately, this assumption is rarely valid for waters with significant particulate (e.g., algal and mineral) backscattering [2]. If R_{rs}(NIR) can be modeled, however, it can be removed from the TOA signal, allowing the atmospheric correction to proceed in estimating ρ_{a}(NIR) as if R_{rs}(NIR) were negligible.

The OBPG implemented a bio-optical model to estimate R_{rs}(NIR) as part of the third reprocessing of the SeaWiFS mission in May 2000 (R2000) [2,3]. This approach required an initial run through the atmospheric correction process (ignoring R_{rs}(NIR) > 0 sr^{−1}) to estimate C_{a}. This preliminary C_{a} was used to estimate spectral particulate backscattering, b_{bp}(NIR), which was used in turn to reconstruct R_{rs}(NIR). With this modeled R_{rs}(NIR) removed from the TOA signal, the atmospheric correction process was repeated, C_{a} was recalculated, the process was iterated upon until convergence in C_{a} was achieved [2,4]. With the fourth SeaWiFS reprocessing in July 2002 (R2002), this C_{a}-driven algorithm was replaced with a reflectance-based algorithm [5]. The OBPG switched from [2] to [5] after observing that [2] depressed ratios of ρ_{a}(NIR) in optically complex waters. This often resulted in the selection of spectrally flat aerosol models within the atmospheric correction process, which artificially depressed the final C_{a} retrievals. The initial estimates of C_{a} were unrealistically high in these complex water masses, which led to over-estimates of b_{bp}(NIR) and, subsequently, R_{rs}(NIR). The reflectance-based approach [5] was less prone to such over-correction. In this revised model, R_{rs}(λ) from an initial run through the atmospheric correction process and empirical estimates of the absorption coefficients for particles and dissolved materials at 670-nm were used to estimate b_{bp}(670) via an ocean surface reflectance model [6]. A spectral scattering function [7] was then used to derive b_{bp}(NIR) from b_{bp}(670). Finally, R_{rs}(NIR) was reconstructed from b_{bp}(NIR), the atmospheric correction process was repeated, and the process was iterated upon until convergence in R_{rs}(NIR) was achieved.

Despite such efforts to account for non-negligible R_{rs}(NIR), the atmospheric correction algorithm continued to demonstrate sub-optimal performance in highly scattering waters, as manifest by negative R_{rs}(λ), particularly in the blue region of the spectrum, and erroneous spectral dependence in visible (VIS) part of the spectrum, with the latter leading to inaccurate C_{a} retrievals. To address these issues, a re-evaluation of the implementation of [5] was initiated prior to the sixth SeaWiFS reprocessing in September 2009 (R2009) and several modifications were made to the R_{rs}(NIR) model, which are presented below. First, an inconsistency in the relationship between R_{rs}(λ) and the inherent optical properties (IOP) was eliminated. Second, alternative estimates of the spectral dependence of b_{bp}(λ) and of the absorption coefficients for particles and dissolved materials at 670 nm were adopted. Finally, a revised iteration scheme was implemented.

## 2. NIR reflectance model development

**2.1** Inversion of remote-sensing reflectance

Remote sensing reflectance is proportional to the IOPs, specifically the marine absorption and backscattering coefficients, through:

whereG(λ) is a function of the illumination conditions, water constituents, and sea-state, b_{b}(λ) is the total backscattering coefficient [ b_{bw}(λ) + b_{bp}(λ)], and a(λ) is the total absorption coefficient [ a_{w}(λ) + a_{p}(λ) + a_{g}(λ)]. The additional subscripts w, p, and g indicate specific contributions by water, particles, and dissolved material, respectively.

The NIR reflectance model [5] implemented for R2002 and used through SeaWiFS reprocessing 5.2 in 2007 (R2007) suffered from an inconsistent application of Eq. (1). First, an approximation of this relationship, where G(λ) was fixed at 0.51 and X(λ) was simplified to b_{b}(λ) / a(λ), was used to obtain near-infrared IOPs from R_{rs}(VIS). The quadratic form of [6], where G(λ) = [g_{1} + g_{2} X(λ)], g_{1} = 0.0949, and g_{2} = 0.0794, was later used to obtain R_{rs}(NIR) from the modeled near-infrared IOPs. While a minor inconsistency, using dissimilar definitions for G(λ) and X(λ) was found to have significant effects in highly scattering waters.

Furthermore, the approach was inconsistent with the methods of [9] that are used elsewhere (and independently) in the atmospheric correction process to normalize R_{rs}(λ) for bi-directional effects. In [9], G(λ) = *f*(λ) / Q(λ), where *f*(λ) is a function of the illumination conditions, water constituents, and sea-state, Q(λ) is the ratio of spectral upwelling irradiance to spectral upwelling radiance, and ratios of *f*(λ) to Q(λ) are pre-tabulated from extensive radiative transfer simulations. For the R_{rs}(NIR) model implemented in R2009, G(λ) = *f*(λ) / Q(λ), using the look up tables (LUTs) of [9], for both the conversion of retrieved R_{rs}(VIS) to IOPs and the conversion of modeled IOPs to R_{rs}(NIR). While these LUTs were originally generated for waters whose light field is predominantly influenced by phytoplankton, they proved suitable for use in this model, as will be demonstrated in Section 4. A full discussion of the applicability of [9] to all water types is beyond the scope of this paper.

**2.2** Spectral backscattering

Given the relationship in Eq. (1), a measured R_{rs}(λ_{0}), and knowledge of a(λ_{0}), an estimate of b_{b}(λ_{0}) can be obtained, where λ_{0} is any reference wavelength (e.g., 670 nm). This retrieved b_{b}(λ_{0}) can then be extrapolated to the NIR through a spectral model. The spectral dependence of b_{b}(λ) can be expressed as:

In Eq. (2)b, the η term is typically in the range of [0, 2], with values around 2 observed in open ocean waters and turbid waters approaching zero [12]. In contrast, the model employed in R2002-R2007 [3] used a spectral scattering relationship for Eq. (2)a developed by [7]:

_{rs}(NIR) model, the latter relationship required an assumption of equivalent spectral dependency for backscattering and scattering. In the R

_{rs}(NIR) model employed in R2009, the more common power-law form of Eq. (2)b was adopted. Rather than relying on a fixed η value, however, the empirical relationship developed in [12] was adopted:

The NASA bio-Optical Marine Algorithm Data set (NOMAD) [13] includes measurements of b_{b}(VIS) at a number of wavelengths spanning a range of values from 0.0007 to 0.0133 m^{−1}. Unfortunately, this data set does not include measured b_{b}(NIR). Nevertheless, this data set was used to evaluate the various forms of Eq. (2) and 3 using a reference wavelength of 670 nm and a retrieval wavelength of 555 nm. The formulation of [12] had the minimal bias and median percent difference relative to in situ measurements (Table 1
, Fig. 1
). and performed better than the linear, spectral scattering function [7].

**2.3** Estimating absorption

In the NIR, a(λ) is dominated by the absorption due to water, a_{w}(λ) [10,11], such that absorption coefficients for particles, a_{p}(NIR), and dissolved material, a_{g}(NIR), can be ignored. However, a reasonable estimation of the absorption coefficient is crucial for the inversion of R_{rs}(λ) to retrieve b_{b}(670), where a_{p}(670) and a_{g}(670) can be a significant portion of the total. In [5], the estimate of a(670) involved two functions, one based on C_{a} [14] and the other on a green-to-red reflectance ratio [5]. In the modified R_{rs}(NIR) model, this was updated to a single, C_{a}-based relationship derived from NOMAD [13] (Type II linear regression using spectroscopy only, limited to C_{a} > 0.2 mg m^{−3}; Fig. 2
):

_{p}(670) and a

_{g}(670) values contained within the NOMAD data set span a range of 0.00001–0.03 m

^{−1}for a

_{p}(670) and 0.00057 0.8413 m

^{−1}for a

_{g}(670).

## 3. Iteration scheme

Since the retrieval of R_{rs}(λ) in highly scattering waters requires knowledge of R_{rs}(NIR), and estimation of R_{rs}(NIR) requires knowledge of R_{rs}(VIS), the R_{rs}(NIR) algorithm must be implemented iteratively. The iteration begins by assuming R_{rs}(NIR) = 0 sr^{−1}, so that initial R_{rs}(λ) can be derived and a first estimate of R_{rs}(NIR) can be achieved. The iteration converges when R_{rs}(λ) changes by less than 2%. In the majority of cases, the method converges after 3-4 iterations, but the algorithm allows up to 10 iterations.

If, however, the initial atmospheric correction results in negative or an otherwise non-physical spectral distribution for R_{rs}(λ) (criteria based on in situ observations), the iteration is re-initialized assuming ρ_{a}(NIR) = 0 (i.e., all reflectance except that from Rayleigh scattering and Sun glint originates from the water mass). This is effectively the opposite extreme from the initial condition of R_{rs}(NIR) = 0 sr^{−1}. The iteration process is then allowed to proceed to convergence as previously described. For the rare cases where the iteration fails to converge after this reset, an additional iteration is forced, again with an assumption of ρ_{a}(NIR) = 0. This pixel is then flagged with an ‘atmospheric correction warning’ to alert users to the questionable nature of the pixel, although a qualitatively useful retrieval may still result.

The controls on this iteration scheme differ slightly from previous implementations [8], especially with regards to the how the iteration is reset, which is currently more robust than that of past methods. It should also be noted that the R_{rs}(NIR) is forced to zero if the initial iteration results in a C_{a} < 0.3 mg m^{−3}. This is done to avoid the addition of noise into clear-water pixels through the estimation of otherwise negligible R_{rs}(NIR) contributions. Furthermore, to avoid the introduction of transitional artifacts, the R_{rs}(NIR) estimation is linearly weighted from 0 to 1 for 0.3 ≤ C_{a} ≤ 0.7 mg m^{−3} (Fig. 3
). This C_{a}-based weighting of R_{rs}(NIR) was first introduced in R2002 and is effectively unchanged for R2009.

## 4. Impact

A long-term satellite and in situ data set was recently compiled to study C_{a} variability in Chesapeake Bay [15]. These data were used to evaluate the modifications to the NIR correction described above, as implemented for the most recent SeaWiFS reprocessing. Specifically, the NIR correction implemented with R2002 (and used through R2007) was compared to the method described in this work and employed for R2009. The satellite and in situ data acquisition and processing methods are described in [15].

The revised correction reduced the percentage of negative R_{rs}(412) for the time series of SeaWiFS retrievals over the Bay by nearly half (18.5% vs. 31.9%). The incidence of negative R_{rs}(490) was reduced to negligible levels (Table 2
). The satellite estimates of C_{a} showed a substantial improvement in the agreement with ground truth measurements in the turbid and highly productive regions of the Middle and Upper Bay (Fig. 4
).The revised correction also improved reflectance retrievals relative to coincident in situ measurements. Following the method of [16], and using a global validation data set compiled by the OBPG, we compared R_{rs}(λ) retrievals for SeaWiFS after limiting the data set to where the NIR correction was used in the processing of the satellite data (e.g., C_{a} ≥ 0.3 mg m^{−3}). Improvement in the satellite agreement with in situ R_{rs} was evident for all bands, with the greatest improvement seen for the 412 and 443 nm bands (Table 3
).

## 5. Conclusion

The revised R_{rs}(NIR) model that is presented in this paper has been implemented in the operational processing of satellite ocean color sensor data by the NASA Ocean Biology Processing Group. This revised correction eliminated a number of inconsistencies present in the previous version. A variable estimate for the spectral slope of the backscattering coefficient was implemented in place of a fixed linear relationship. A simplified estimate of absorption in the red region of the visible spectrum was also implemented. This resulted in a significant reduction in the incidence of negative R_{rs} for satellite retrievals in the 412 to 490 nm range, while improving the agreement of satellite-derived R_{rs}(λ) and C_{a} with in situ measurements.

The results and discussion presented here have been focused on ocean color retrievals from the SeaWiFS sensor. It should be noted, however, that the same algorithm has been applied in the processing of MODIS-Aqua and other ocean color sensor data sets archived by NASA’s Ocean Biology Processing Group.

The approach presented, while a significant improvement over the previous implementation, cannot correct all cases where the black pixel assumption is invalid. The use of empirical optical models (e.g. spectral backscattering, particulate and dissolved absorption) limits the applicability of this algorithm to waters that are similar to those over which these empirical models were developed. In addition, the practical implementation of the iteration scheme can result in conditions of non-convergence. Fortunately, the optical domains covered by the data sets used in the generation of these algorithms are quite diverse and the cases of convergence failure are small. Therefore, these limitations do not detract from the benefits realized by the use of this algorithm.

## References and Links

**1. **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] [PubMed]

**2. **D. A. Siegel, M. Wang, S. Maritorena, and W. Robinson, “Atmospheric correction of satellite ocean color imagery: the black pixel assumption,” Appl. Opt. **39**(21), 3582–3591 (2000). [CrossRef]

**3. **C. R. McClain, E. Ainsworth, R. Barnes, R. Eplee, Jr., F. Patt, W. Robinson, M. Wang, and S. W. Bailey, “SeaWiFS Postlaunch Calibration and Validation Analyses, Part 1,” NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2000).

**4. **A. Morel, “Optical modeling of the upper ocean in relation to its biogenous matter content (case 1 waters),” J. Geophys. Res. **93**(C9), 10749–10768 (1988). [CrossRef]

**5. **R. P. Stumpf, R. A. Arnone, J. R. W. Gould, P. M. Martinolich, and V. Ransibrahmanakul, “A partially coupled ocean-atmosphere model for retrieval of water-leaving radiance from SeaWiFS in coastal waters,” in Patt, F.S., et al., 2003: *Algorithm Updates for the Fourth SeaWiFS Data Reprocessing*. NASA Tech. Memo. 206892, National Aeronautics and Space Administration, Goddard Space Flight Center, Greenbelt, MD (2003).

**6. **H. R. Gordon, O. B. Brown, R. H. Evans, J. W. Brown, R. C. Smith, K. S. Baker, and D. K. Clark, “A semianalytic radiance model of ocean color,” J. Geophys. Res. **93**(D9), 10909–10924 (1988). [CrossRef]

**7. **R. W. Gould Jr, R. A. Arnone, and P. M. Martinolich, “Spectral dependence of the scattering coefficient in case 1 and case 2 waters,” Appl. Opt. **38**(12), 2377–2383 (1999). [CrossRef]

**8. **NIR, html, NIR Correction, http://oceancolor.gsfc.nasa.gov/REPROCESSING/SeaWiFS/R4/NIR.html

**9. **A. Morel, D. Antoine, and B. Gentili, “Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function,” Appl. Opt. **41**(30), 6289–6306 (2002). [CrossRef] [PubMed]

**10. **L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65-2.5 m spectral range,” Appl. Opt. **32**(19), 3531–3540 (1993). [CrossRef] [PubMed]

**11. **R. M. Pope and E. S. Fry, “Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,” Appl. Opt. **36**(33), 8710–8723 (1997). [CrossRef]

**12. **Z. Lee, R. A. Arnone, C. Hu, P. J. Werdell, and B. Lubac, “Uncertainties of optical parameters and their propagations in an analytical ocean color inversion algorithm,” Appl. Opt. **49**(3), 369–381 (2010). [CrossRef] [PubMed]

**13. **P. J. Werdell and S. W. Bailey, “An improved in-situ bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. **98**(1), 122–140 (2005). [CrossRef]

**14. **A. Bricaud, A. Morel, M. Babin, K. Allali, and H. Claustre, “Variations of light absorption by suspended particles with the chlorophyll a concentration in oceanic (case 1) waters: analysis and implications for bio-optical models,” J. Geophys. Res. **103**(C13), 31033–31044 (1998). [CrossRef]

**15. **P. J. Werdell, S. W. Bailey, B. A. Franz, L. W. Harding Jr, G. C. Feldman, and C. R. McClain, “Regional and seasonal variability of chlorophyll-a in Chesapeake Bay as observed by SeaWiFS and MODIS- Aqua,” Remote Sens. Environ. **113**(6), 1319–1330 (2009). [CrossRef]

**16. **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]